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We show that rhombohedral graphite may support surface superconductivity with an unusual 
relation between the BCS coupling constant and the order parameter. This feature results from 
the properties of the states localized on the graphite surfaces. In a description including only the 
nearest neighbour coupling of the graphene layers, the surface states are topologically protected 
and have a flat band dispersion. We show that including higher order couplings destroys this flat 
band character and leads to a particle-hole symmetry breaking quadratic dispersion with a large 
effective mass. Employing this dispersion, we then show its effect on superconductivity and find two 
regimes of parameters, depending on the relation between the strength of the coupling constant and 
the details of the quadratic dispersion. For low coupling strengths, superconductivity is localized 
on the surfaces, but the order parameter is exponentially suppressed as in a conventional BCS 
superconductor, whereas for large coupling strengths we obtain surface superconductivity with a 
linear relation between the order parameter and the coupling constant. Our results may explain the 
recent findings of graphite superconductivity with a relatively high transition temperature. 

I. INTRODUCTION 

Superconductivity is a ubiquitous phenomenon in metals: according to a commonly held view, all metals become 
either superconducting or magnetic at low enough temperatures. However, the corresponding transition temperatures 
may be (so far) unobservably low. In conventional superconductors, such as Al, Hg, or Nb, the transition temperature 
depends exponentially on the inverse of the BCS coupling constant, and whereas the coupling constant itself may be 
fairly large (the relevant energy scale connected with it may be many times larger than the thermal energy at room 
temperature), the resulting transition temperature typically does not much exceed 10 K. This property is intrinsically 
related to the quadratic dependence of electrons' energy on momentum, which leads to a logarithmic divergence in 
the BCS self-consistency equation. With a higher-order dispersion around the Fermi energy, the relation between the 
magnitudes of the critical temperature and the coupling constant becomes stronger, boosting the superconductivity. 

The extreme case would be a completely dispersionless energy spectrum, the so called "flat band". Fermionic 
systems with dispersionless branches of excitation spectrum have quite unusual properties; nowadays they attract lots 
of research interest. Flat bands were predicted in many condensed matter systems, see for example^^. In some cases 
the fiat bands are protected by topology in momentum space; they emerge on the surfaces of gapless topological matter^ 
such as surfaces of nodal superconductors"^, graphene edges', surfaces of mult ilayere d graphene structures'^, and 
in the cores of quantized vortices in topological superfluids and superconductors^ 13 1 14 1 . The singular density of states 
(DOS) associated with the dispersionless spectrum was recently shown by us to essentially enhance the transition 
temperature opening a new route to room-temperature superconductivity. 

The problem is to f ind th e metal with such a higher-order dispersion around the Fermi sea. Along with our collabo- 
rators, we have showrP^E^ that within the nearest-neighbour approximation, rhombohedral graphite has topologically 
protected surface states with a flat band at the Fermi energy, and these surface states support high-temperature 
superconductivity where the superconducting order parameter is concentrated around the surfaces. Such a supercon- 
ductor may also carry a large surface supercurrent with a critical value proportional to the large critical temperature. 
The corresponding critical temperature depends linearly on the pairing interaction strength and can be thus consid- 
erably higher than the usual exponentially small critical temperature in the bulk. A flat band forms out of a low 
dispersive band that appears on the surface of a multilayered graphene structure with rhombohedral stacking with a 
large number of layers. Surface superconductivity is favorable already for a system having TV > 3 layers, where the 
normal-state spectrum has a power-law dispersion £ p oc |p| as a function of the in-plane momentum p. The DOS 

j/(£ p ) oc £p 2 N " N has a singularity at zero energy which results in a drastic enhancement of the critical temperature. 

However, next-nearest neighbour hoppings which are present in real rhombohedral graphite can break the exact 
topological protection and, therefore, the flat-band mechanism of superconductivity at sufficiently low values of the 
coupling constant can be destroyed. Here we study the detailed effect of these higher-order interactions and show that 
though they indeed break the flat-band scenario for weak superconducting coupling, they provide another mechanism 
of surface superconductivity which is of the BCS type but has a much larger coupling constant than the usual 
superconductivity in bulk graphite. This large coupling constant comes from a large DOS associated with a heavy 
effective mass of surface quasiparticles that is clearly distinguishable on the background of the flat band which would 
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exist without the higher-order interactions. Both these mechanisms favor the high-temperature superconductivity. 
Our results provide a criterion for the parameters needed to obtain the highest critical temperature. They may 
be relevant in explaining the recent experimental findings^"^ reporting the observation of even room-temperature 
superconductivity in doped graphite. 

Here we briefly outline the main results of our paper. More detailed calculations are described in the following 
sections. In this work we study rhombohedral graphite by taking into account, besides the lowest-order interlayer 
hopping energy 71, also the higher-order hoppings 73 and 74. From these, only the latter breaks the topological 
protection of the flat band. We find that, even in the presence of all these interactions, rhombohedral graphite still 
has surface states. However, instead of a flat band of radius pfb = 7i/Vf they have a weak dispersion 

e p = a[—\ -fi, (1) 
\PfbJ 

within the region p < pfb- Here a small factor a = 27174/70 arises from the higher-order interlayer hopping 74; 
70 is the zero-order intralayer nearest-neighbour hopping energy, vf is the Fermi velocity in graphene, and fi is the 
chemical potential. These surface states are not symmetric with respect to the point e = 0, more or less in the same 
way as the presence of the next-neighbour coupling breaks the electron- hole symmetry in graphene 21 . For p > Pfb, 
the energy of the surface states deviates rapidly from [i (see Figs. [4] and [5]). Qualitatively, the energy gap is localized 
at the surface and can be determined by a simplified self-consistency equation, which for T = and [i = has the 
form 

1 = -\- f " pdp * = /-Arsinh(a/A) , (2) 

where g — Wp FB /h 2 d is a superconducting coupling energy proportional to the pairing interaction W, d is the distance 
between the layers. Equation ^ yields 

A = asinh^ 1 (Aa-rr/g) . (3) 

We thus find that for g < g c ~ 47r7 1 7 4 /7 , the gap A is exponentially small, similar to conventional BCS supercon- 
ductors, whereas for larger coupling constant A tends to the flat-band resullP^ A oc g which is linear in the coupling 
strength. These values of A yield critical temperatures T c of the order A/fee , the exact prefactor depending on which 
of the above regimes the coupling constant is. 

The above considerations help to identify the regime of parameters where extremely high-temperature surface 
superconductivity might be found. In our estimates, we use the tight-binding parameters summarized irP2 (see 
alscP^l), which gives 70 = 3.2 eV, 71 = 0.39 eV, 73 = 0.315 eV, and 74 = 0.044 eV. We disregard the next-nearest 
neighbour intralayer hopping proportional to 72 = —0.02 eV and even smaller interaction across two layers. In 
addition of its small relative magnitude, the next-nearest hopping 72 does not affect the conical dispersion of single 
graphene layer and thus is not expected to essentially modify our results. Note, however, that also widely different 
values of higher-order hopping constants are discussed in literature^, especially with a much higher 74 (for example^ 
claims 74 = 73). With the above values, the crossover between the two regimes takes place around g c w 0.1 . . . O.271 ~ 
0.04. . . 0.08 eV (see Fig. [9] below). Thus, for g < g c , the gap and the critical temperature would be exponentially 
small whereas for g > g c it would be T c ~ (3/71) x 60 K. This is much larger than the expected gap in the bulk for 
the same magnitude of the coupling. 

This paper is organized as follows. In Sec. [II] we derive the dispersion of the surface states for rhombohedral graphite 
in the normal state by including the interlayer couplings 71, 73 and 74. In Sec.lTTllwe derive the Bogoliubov-deGennes 



(BdG) equations for the superconducting state, while Sec. IV describes the surface superconductivity in rhombohedral 



graphite. Besides analyzing the crossover between the exponentially damped superconductivity and strong flat-band 
superconductivity, we review our earlier results on flat band superconductivity in the case of finite number of layers 
and on supercurrent carried by the surface states. In addition to this, we analyze the role of fluctuations around our 
mean-field solution and detail our previous prediction^ of interface superconductivity at twinning layers of graphite 
by showing a few example cases. The paper is concluded by a short summary of our work and a comparison to the 
recent experiments on superconductivity in graphite. 



II. ELECTRON DISPERSION IN RHOMBOHEDRAL GRAPHITE 



The rhombohedral graphite lattice and the tight-binding couplings are depicted in Fig. [T] We denote the layers 
(starting from the bottom) by index n, the position of A atoms inside layer n by Ri^ n , the vectors from the A atoms 
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FIG. 1: (Color online) Rhombohedral graphite lattice, tight-binding parameters and the lattice vectors. 



to the nearest B atoms by Sj (j = 1, 2, 3), and the vector between the layers by d. The latter is strictly in the vertical 
direction, i.e., connects a B atom from layer n to an A atom in layer n + 1 (this is the convention leading to er_ 
coupling on the upper diagonal). 

For this lattice structure wc can write the Hamiltonian 



1=0 



(4) 
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where 

N 

H {0) = -70 EE E [^ t (Ri 1 n)V»(Ri,n + *i)+h.c] 
n=l Ri,„ .3=1,2,3 

JV-1 

ijd) = -71 E E [V'f t (Ri ! n + 5i)^ +1 (Ri,„ + d + 5 1 ) + h.c.] 

n=l R ; „ 
JV-1 3 

= -73 E E E ^( R ..»)t B + i( R ..» + d + h - c -] 

n=l R ijn 3=1 
N-l 3 

ff(4) = -74 E E E K f (Ri,n)^ + l(Ri,n + d + + h.C.] 
n=l Ri_„ j'=l 

N-l 3 

-T* E E E [t Bt ( R ^ + Wn+i(R-i,« + d + «5i + «j) 

n=l Ri,„ i=l 

+ h.c] . 

The sum over R^ n goes o ver all A-atoms in layer n. For the hopping constants we use the most recent data in 
graphite according td 21 | 22 f 70 = 3.2 eV, 71 = 0.39 eV, 73 = 0.315 eV and 74 = 0.044 eV and neglect the intra-layer 
next-nearest neighbour hopping term H^ 2 ' proportional to 72 = —0.02 eV, as well as the even smaller hoppings across 
two layers. 

The unit cell vectors are, see Fig. [2] 

a i = y ( 3 > V^) , a 2 = y (3 , -V$) , a 3 = d + Si, 

where clq is the distance between two carbon atoms while d is the interlayer distance in the z direction. The nearest 
neighbor in-plane vectors are 

*i = a„(l ,0) , 5 2 = y(-l , Vd) , 5 3 = y(-l , . (5) 

The single-layer graphene has Dirac points in the Brillouin zone at 

K=-^(v / 3,l),K' = -^(V3,-l). 
oyoao dv3a 
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In the vicinity of these points the graphene spectrum has a form of touching cones (for details, see revie'wPS and 
references therein). We show below that, for surface states in multilayered graphene with rhombohedral stacking as 
shown in Fig. [T] the conical spectrum near the Dirac points is transformed into low-dispersion low-energy bands that 
determine the unique features of this system. 

Therefore, we are interested here in low energies and thus in momenta close to one of the two non-equivalent Dirac 
corners K or K' in the Brillouin zone. Let us expand the wave functions near K , 

^(RO = 4rE e * (K+p/;i) ' R ^(p) (e) 

v J-/ 

p 

<£(R, +S) = ^=Y1 e^^ +s ^(p) , (7) 
p 

and assume that |p|, |p| <C ha^ 1 . Here L is the number of unit cells in the plane (we omit the spin index). Small 
vectors p cannot couple the two Dirac points K and K'. Therefore, equations for them separate. 
A standard Fourier series expansion of Hamiltonian Q near K yields 

N 

H K = ]T ^(p)4nn(K,p)4(p) , (8) 

p m,n— 1 



H mn (K,p) = 5^ffW„(K,p) , (9) 

where 



2=0 



ffi°i(K,p) = v F (&-p)5 mn 

tf«(K,p) = - 7l [e-™/ 6 a + 6 m , n+ i + e i7r / 6 a_5 m , n _i 

H£l(K,p) = ^v F fe- i7r / 3 (a + p+)5 m , n _ 1 
7o L 

+ e i7r / 3 (a_p_)£ ro , n+1 

fi£l{K,p) = ^« F L" /6 ! ,J ™ J n-i + e" l/ V« m .»+i 
7o L 



and Vp = 3ao 7 o/2S. We define the pseudo-spinor 

& = ( $f) , & = W f , , 

where tp^ = ip^, -0n = e l7T ^ 6 ^^ . Above we use 

2a± =a x ± i&y , p± = p x ± ip y = pe ±l4, . 
We also need the Hamiltonian expansion near the opposite Dirac point — K, 

#*(R0 = -L£ e *C-K+p/»)-R^A (p) (10 ) 
VjL p 

^(R i + 5) = -^^eWP/ft)-^^)^^). (11) 
Using the same type of derivations we find 



L 

p 



AT 

#-k = Yj E ^„(P)^1(-K,p)^(p) , (12) 

Z,p m,n— 1 



where £Z^(-K,p) = H^l*(K,p) and F^(-K,p) = (K,p) for I ^ 1. The Hamiltonian Eq. |l2|) at 

— K together with the Hamiltonian Eq. ^ at K are used below to construct the associated Bogoliubov-de Gennes 
Hamiltonian to describe the superconducting state. 

In the next section we find numerically the energies and eigenstates of Hamiltonian Eq. ^ using the relative 
magnitudes of the coupling constants listed above. The result of such numerics is displayed in Figs. [4] and [5j along 
with the corresponding analytical approximations. 
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A. Low-energy spectrum in the normal state 

The Schrodinger equation takes the form 

^H nm (K,p)^ m (p) = (e + /i)^„(p) . 

m 

The energy E = /i + e is measured from the chemical potential fi. We use the Ansatz 

A 1 



(13) 



where q is the out-of-plane momentum, and obtain for n 7^ 1, N 



v F pe- 1 ^ - lie - lqd + ^v FP e l ^ +qd) 
7o 



74 

2 — vfP cos((p — qd) — (e + /i) 
7o 



7o 



74 

2 — vfPCOs(4> — qd) — (e + //) 
7o 



A 

A 1 =0 
A 1 







(14) 



with <p = <p + 7r/6. The compatibility condition yields the energy spectrum for excitations in the bulk, 

i 2 



74 

e + /i — 2 — vppcos((j) — qd) 
7o 



+ ( — ) M 2 + 2 — M 2 cos(20 + gd) - 2 7l — (vp) cos(</> + 2 9 d) . 



7o 



,73 



u 2 p 2 + 71 — 2(up)7i cos(0 — gd) 



73 



7o 



7o 



For zero doping fi = 0, the Fermi surface is determined by e(p, q, <p) = 0. If 73 = 74 = 0, the Fermi surface shrinks 
to the spiral VfP = 71 , 4> = qd. If only 74 = 0, the Fermi surface for zero doping is determined by 



v F p _ e^-iV + (7 3 / 7 o)e =Fl( * +2 ' 5fd) 
7i 1 + (73/70) 2 + 2(73/70) cos(20 + qd) 

Real momenta p and q are realized along the line q(4>) satisfying 



(15) 




VFPyh 
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FIG. 3: Fermi line of rhombohedral graphite for 74 = 0, 73 7^ (red) compared to the simple spiral obtained for 74 = 73 = 
(blue). 



sin(^ - qd) - (73/70) sin(0 + 2qd) = 



7 



Along this line 



tan 7 = sm(qd) + (73/70) sm(2qd) 
cos(qd) - (73/70) cos(2qd) 



cos(dq — (f>) + 73 cos(2qd + ^ 
7i 1 + 7| + 273 cos(<?d + 20) 



It is a corrugated spiral in the 3D space (p, q, (f>). Since <f> at the Fermi line as a function of qd is periodic with period 
2ir, the corrugation has a three-fold symmetry in the plane of 4>, as is seen in Fig. [3] and in Fig. [4] below for the 
corresponding flat band. It is clear that the 73 interaction does not destroy the Fermi line. Indeed, this is because 
the interaction comes with matrices a x and a y and thus the full Hamiltonian obeys the same anti-commutation rule 



a z ,(H (0) +£(3)) 



= 



as the initial Hamiltonian + H^> . According tcP^, this preserves the same topological invariant and hence the 
topology of the Fermi surface is unchanged. 




vfPxHi 



FIG. 4: Spectrum of rhombohedral graphite surface state as a function of momenta p x and p y . On this scale, the effect of 74 
does not show up. Due to the non-zero 73, the shape of the low-energy region is distorted from a circle into a more triangular 
shape. 



We now turn to the surface states with low energies e = 0. A surface state corresponds to a complex out-of-plane 
momentum q — q' + iq" that ensures its decay into the bulk. Since the interaction proportional to 73 does not change 
the conclusion about the presence of a flat band, we restrict our analytical consideration to the case when only 74 is 
nonzero while 73 = for simplicity. Note that our numerical analysis is carried out using the exact diagonalization of 
the full Hamiltonian (with non-zero 73). For 73 = Eq. (13) takes the form 



v F (a ■ p)V>„(p) - 71 



74 

7o 



t/6 



(16) 



Equation ( 16 ) suggests a solution in the form 



= e i(0- w / 6 )(n-l-f ) 



A 4 



N- 



i<l> 



(17) 



where 



np \ 71 (e + M) - (74/70) («V + ll\ 
7i, 



v 2p2 



7? 




FIG. 5: Cuts of the 3d spectrum shown in Fig. [4] along the p x (solid lines) and p y (dash-dotted lines) directions and limited 
to the low-energy region vfp < 71. Three cases are shown, TV = 5 (blue), TV = 20 (red), and TV = 50 (black) graphene layers. 
The inset shows a zoom-up of the low-energy region. The dashed lines show the corresponding approximations from Eq. ( 20 1 , 
without including the correction l — v 2 p 2 /j 2 as this becomes relevant only in the superconducting state. The deviations between 
the dashed and other lines show up mostly when the term £ p becomes dominant in Eq. (20 1, and are partially due to the term 
73 neglected in that approximation. 



This form of solution implies the out-of-plane momentum q'd = <j> while e ±q d — (vfp/ji)- 

At the outermost layers n = 1,TV, the terms which would contain tpQ and tpN+i m Eq. (16) disappear. Those 
components for which the terms with 71 are absent have the form 



7o 



74 



7o 
vppe 1 



-icj>+iir/6^pl 



= (e + fijipl 



N-l 



N 



(18) 
(19) 



They couple the constants A + and A and determine the energy of the surface states. We define 

C P = 7i (vfp/ii) N ■ 



Using Eq. ( 17 1 we find for £ p , e -C 71 
The normal-state dispersion is e = e„, 



74 



ll 



2 2 (e + M)-27i 2 

v z p 2 70 7f 



2 2 
v z p z 



v 2 p 2 



A 



fj, p ± C P (l - «V/7i) > 



2m* 



where 



7i 7o 



9p 2 474WJ, 



(20) 



(21) 



The spectrum e p has a quadratic dispersion with the effective mass m* on a background of a much weaker high-order 
dispersion £ p . The latter transforms into a flat band £ p = with a radius vfp < -fi for an infinite number of layers, 
TV — > 00. This form of the dispersion is compatible with the findings made with the numerical calculations exploiting 
density functional theory also there quadratic dispersion with an effective mass of the same order as that estimated 
here is found. 

The effective mass is much larger than the characteristic band mass in 3D graphite. Indeed, we have 



m 
m 3 



71?* 2 



7470771305 



71 

74 
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where we estimate ti 2 /(m^a^) ~ 70 as the conduction band width in graphite. We see that m* /m.3 3> 1. The group 
velocity is v g = de p /dp = p/m* . 

This dispersion is compared with the exact dispersion obtained from the numerical diagonalization of H(K,p) in 
Figs, [i] and [5] using 73/70 = 0.098, and 74/70 = 0.014 according tcP^S 



III. BOGOLIUBOV-DE GENNES EQUATIONS FOR THE SUPERCONDUCTING STATE 



To construct the BdG equations we need the Hamiltonian for holes J"™ The hole Hamiltonian near a Dirac point K 
follows from the particle Hamiltonian in a vicinity of the opposite Dirac point — K. The wave function ip^ of a hole 
excitation near the Dirac point K is V>k = ^-k- Therefore, the hole Hamiltonian is H h (K, p) = H*(— K, — p). Since 



the term Hmh(K-, p) is independent of p while the other terms are linear in p we have using Eq. ( 12 ) 



N 



f k = E E ^ )t (p)^(-K,-p)^ h) (p) 



/,p m,n=l 
N 



E E ^HpMUk^^Hp) 



l,p m,n— 1 
N 



p m,n— 1 



E ^ )t (P)^m„(K,p)^( P ) 



In other words, the hole Hamiltonian coincides with Eq. (|9| where the electronic wave functions are replaced with 
the corresponding hole functions. As distinct from the quasiparticle energy measured from the chemical potential 
upwards, E = /1 + e, the energy of holes is measured from the chemical potential downwards, E = ji — e. In what 
follows, the electron wave function is denoted by u n = ip n while the hole wave function is denoted by v n = tpn- 

The BdG equations for the superconducting state are constructed using the Schrodinger equations of the type of 
Eq. (16) for particles and holes with the particle-hole coupling through the order-parameter field A, 



E f 3 ® H nm (K,p) - fi5 

nr. 



A 



(22) 



Here we introduce objects in the Nambu space 



T3 = 



1 

-1 



A„ = 



A n 
At 



The Nambu vector ty n has the pseudo-spinor components 



For 73 = the BdG equations in components take the form (for n 7^ 1,N) 



v F (a-p)u n (p)- 7l (W 6 a 



-U n +1 + e " /6 (T + U n _i 



74 / i7r/6 - 1 -in/6 

( e 1 v F p-U n+ i + e ' v F p+u n -i 



7o 



//(/ 

+A n v n = eu n , 
v F (& ■ p)v n (p) - 71 (e' i7r/6 <r_{i n+1 + e-^^a+in-! 



^ (e™/ 6 v F p-i>n+i + e-"/ 6 v FP+ v n ^ 
7o v 



(23) 



(24) 



At the outermost layers n=l,N, the terms with uq, vq and wjv+ij ^iv+i in Eqs. ( 23 ) , (24| disappear. For those 
components which do not contain the terms with 7! we have 

74,. _-i0+ 47 r/6Q,+ 



v F pe ^T^d-, + —v F pe 
7o 



— fia^ + Aia^ = 



v F pe v r 3 a^ 



74 

— \ 
7o 



-v F pe 



i(j>-iir/6-- 



& n-i ~ M"at + A N a 



N 



CO' 



N 



(25) 
(26) 
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Here we decompose the wave function 



* = 



ft 



® V" 



(27) 



into the spinor functions localized at each sublattice 



and introduce the vector in the Nambu space 



For analytical consideration we assume that A„ = for n =/= 1,N. This assumption is justified by the results 
of numerical solution of the self-consistency equation with the wave functions found from the full BdG equations, 
discussed below, see especially Fig. 10 In this case, Eqs. ( 23 1— ( 24 ) for n ^ 1,N do not contain A, so that one can 
use the normal-state coefficients as in Eq. (17), 

C 



a, 



C 

7T 



,<(n-l-f)(0-f) 



i(n-l-f )(0-f ) 



(if)' 
(if) 



A+ 



N-n 



A~ 



( vfP 
\ 7i 

if) 



N- 



n-l 



a" 



Here C is a normalization constant; the vectors A ± = (A ± 1 B ± ) T do not depend on n. We also define the matrix 

'vp\ 7i (f 3 e + n) - (74/70) {v 2 p 2 +7i) 



C 



7i 



^2^2 _ ^,2 



Equations (25 )-(26 1 yield 



where 



f^ p A = (e - f 3 /2 p )i+ - Ail+ , 
-T3? P i + = (e - r 3 /i p )i _ - AjvA- , 



(28) 
(29) 



e = e(l-uV/7i) 1 , = (1 - «V/7i 2 ) 1 • 
Equations (28), (29) provide the surface-state spectrum^ 

[2e 2 + 2^ 2 _ At A N - A 1 A* N ] = . 

If Ai = A^v we have 

~e 2 = (f, p ±Q 2 + \A\ 2 . 

This is compared to the exact diagonalization of the Bogoliubov-de Gennes equations in Fig.JHj 

Equations (28), (29) determine four independent states. If 74 = [/, = they are: (i) e% — E and Af = u, Bf = v, 
(ii) e 2 = — E and .Af = v > ^2 = ( m ) ^3 = & an ^ A* = -B* = ± u , (i y ) ?4 = ~^ an d A^ = ±u, Bf — qpu. 
Here £7 = W^ 2 + A 2 and 



(30) 
(31) 



V2 



1 + 



1 

75 



1 - e P /£ 



The overall normalization requires dJ2n=i[\ a n | 2 + \Pn P + l a « | 2 + IAi | 2 ] = 1- F° r Cp 7i this gives 

iq^d-Hi-KWTi) 2 ] . 



(32) 



(33) 
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-0.5 0.5 

VFPx/y/jl 

FIG. 6: Spectrum in the superconducting state, assuming a coupling constant g = 7i/2 yielding a gap Ai = Ajv = O.OI571 for 
TV = 20 layers. The Bogoliubov-de Gennes equation has four branches of solutions, two for the electrons and two for the holes. 
The direction p y = is shown with the solid lines and p x — with the dash-dotted lines. The black lines are the corresponding 
approximations from Eq. (31 1 and the deviations with the exact numerics are mostly due to the 73 term neglected in Eq. (31 1. 



If the number of layers N is large, £ p — > for vfP < 7, the two surface states decouple 

if = A' + |A|? , e 2 N = /ip + . 



In this case, Eq. ( 28 ) yields 



(e-f 3 A P )i + - Aii + = 
whence A + — U, B + = V or A + = V, B + = U where 

u = ^=[i + n P /iY 1 v = ~^ [1 — h/t 2 ■ 

In the following section we use these eigenfunctions to calculate the self-consistent gap function A. 

IV. SURFACE SUPERCONDUCTIVITY 



(34) 
(35) 

(36) 



The surface states discussed above form a basis for the superconducting order parameter localized near outer 
surfaces. Within the mean-field approximation, the order parameter at layer n is determined by the self-consistency 
equation, 



A„ = 



d 2 p 
(2tt/i) ; 



Y,WTr[u n (p,k)v* n {p,k)\ 



x [1-2/(^,0] , 



(37) 



where W is the 3D coupling potential, f(E) is the Fermi distribution function. The sum includes all states q with 
given 2D momentum p. As was shown irPSI, the terms corresponding to the surface states dominate in the sum due 
to a much larger density of states compared to the bulk states. For a large number of layers N ~ > 00 when £ p = 0, 
the self-consistency equation for the order parameter on the surface takes the form 



A 



-2W I jf^ iqWtanh-1; , 



(38) 
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where pfb = Jx/vp. Using U and V from Eq. (36) 

A = W 



cPp (l-«V/7i 2 ) 2 A 



< l ■I p <pfb ( 2irh ) 2 Jp% + |A| 2 (l-wV/7i 2 ) 2 
v /^ + |AP(l- u V/7i 2 ) 2 



x tanh ■ 



2T 



(39) 



In the following sections we consider several examples that can be derived from Eqs. ( 37 1 — (391 



Flat band 



For zero doping /i 
(for T = 0)P 



0, in the absence of a band curvature 74 = and £ p = 0, Eq. (39) yields the flat-band result 



A ■ 



W 



P<Pfb 



d 2 p 

(2ttH) 2 



1 - 



9 2 



9_ 
8tt 



(40) 



where the coupling energy is 



(W/d)jf _ (W/d)p 



FB 



h 2 



The coupling energy can be expressed in terms of the usual BCS coupling constant A = v$W where u 3 = m^pj,F jlic'T? 
is the 3D density of states. Assuming the conduction band width in 3D graphite of the order of 70 we have 



7? 



3~A 

7o P3Fa 

which can be estimated as 5/71 ~ A (71 /7a) if h/a^p^p ~ 1- 

The critical temperature is determined by Eq. ( 39 ) with A - 



(41) 



0, which gives Aq = 3ksT c . Due to its linear 



dependence on the interaction strength, the critical temperature is proportional to the area of the flat band and can 
be essentially higher than that in the bulk. 
Doping in the flat band regime destroys the surface superconductivity^. This can be seen from Eq. ( 38 ) with 



UV = A/2?i and where ?i is taken from Eq. (34). Both A and T c vanish at the critical doping level \fi\ = 2ksT c 

For a flat band £ p — with p c = pfb the only characteristic values in the superconducting surface state are the 
energy A and the momentum j»fb ■ Therefore, the coherence length should be of the order of the only available length 
scale, £0 ~ fr/PFB- It is much larger than the interatomic distance, £0 3> ao> since pfb po ~ H/ao. 



B. "Flat-band" surface superconductivity in a finite array. 



Let us discuss the "flat band" regime pL 
normal-state DOS defined as 



and 74 = for a system with a finite number of layers. Since the 



p dp _ 7i(f p /7i)' 



2nh 2 d£ p 



2nh 2 Nv 2 F 



(42) 



has a low-energy singularity for N > 2, the surface superconductivity is favorable already for a system with a finite 
number of layers N > 3. A simple expression for the zero-temperature gap can be obtained if iV > 5. For a finite 
N, the value £ p can reach values larger than A. We insert Eqs. (32) at zero doping for u and v together with Eq. 
(|33) and Eq. f3T) into Eq. ((38}. 
C = 71 > A. 



Since the upper limit of integration is Pfb, the corresponding upper limit for £ p is 
Transforming to the energy integral with the normal-state DOS Eq. (42) we find for T = 



1 



W 



^p) 



2Trh 2 dNv 2 F 



d£, p 

^1+ IAI-' 

- (W7.' ^[i-(W7i)^: 

§ + |A| 2 



d£, P 



(43) 
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We see that, for N > 4, the integral converges at £ p ~ A or p ~ p A = p FB (A/7i) « . The zero-temperature gap is 



A = 71 



4"7T7l 



a(JV)- -(Ao/71)" a(iV/2) 



(44) 



where 



a(JV) 



JV + 2 

x ^ ax 
^fx^+f 



N 



2N 



N + 1 

iV 



For iV > 1 we have a^r = 1. The flat-band result, Eq. (40), is recovered if the number of layers is N 21n(7i/Ao). 
The coherence length for a finite system is £0 ~ tl/PA- It approaches H/pfb for N — > 00. 

The gap obtained by numerical integration of Eq. (37 1 with a cut-off p co is plotted in Fig. [7] It shows how the gap 
is independent of the number of layers as long as A > £ Pc0 , and saturates for £ Pe0 > A. As the number of layers N or 
the coupling energy g increases, this threshold moves to larger values of cutoff momenta. 
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FIG. 7: (Color online) Zero-temperature gap as a function of the momentum cutoff p c for various N (solid lines). The gap 
saturates at p c ~ pa and approaches Eq. (401 for TV — > 00. The dashed lines show the dispersion £ p for each TV. Here g = O.OI71 



and the higher-order couplings 73 and 74 have been set to zero. 



C. BCS-like surface superconductivity for a quadratic spectrum 



In this section we consider a system with an infinite number of layers, such that £ p = 0. The spectrum has a weak 
dispersion Eq. (31 1 due to a large effective mass m* determined by the inter-layer coupling 74. In this case Eq. (39) 
yields for T = 



9 

v g 



d 2 P 



7l Jp<p FB (2?0 2 ^2 + ^12(1-^2/^)2 ' 

A simple qualitative expression can be written neglecting the normalization factor 1 — v 2 p 2 /jf: 

d 2 p 1 



(45) 



1 = 



v 2 g 



7? Jp< PFB (2tt) 2 O + | A | 2 



9 
Ana 



Arsinh 



a — ji 



Arsinh- 



14 



where a = 271(74/70). For /j, = and fi = a we find 



A 



a 



smh(Aita/ g) 



This is Eq. ([2| discussed in Sec. [T] In the limit j » a we get the flat band result A = g/4n (without the extra 
1/2 because of the absence of the normalization factor 1 — v 2 p 2 /j±). Whereas this function captures the qualitative 
aspects of the self-consistent gap (exponential suppression below g < Aira, and flat-band limit for g > 47ra), the exact 



numerical value for the gap needs to be found from the full equation ( 45 1 
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FIG. 8: Optimal doping level vs. the coupling constant g for a given value of a. Inset shows the dependence of the (normalized) 
self-consistent gap on the chemical potential for a few values of g. In the limit g 3> 47ra, the scale for the /i-dependence is 
A > a rather than a, but the optimum is still found around /1 « 0.7a. 

For a weak interaction A <C a the gap Eq. ^ has a BCS-like form A = ae -1 ^ 2 where the coupling constant is 
A2 ~ 9/0- Using the estimate Eq. (41) for g we find that this is a much larger coupling constant 



A 2 - A(7i/7 4 ) > A 

than what one would have for the usual bulk superconductivity. The coherence length has its usual form 

Co = HVg/A - Hp FB /m*A - a (7o/7i)e 1/A2 , 

which is much longer than the interatomic distance ao- 

Inserting a non-zero chemical potential fj, in Eq. ( 45 1 enables us to find an optimal doping with which the gap A is 
maximized. Such an optimal doping depends on the ratio g/a, vanishing at g -C a and saturating to a finite value 
fi opt « 0.7a at g 3> a. The optimal doping is plotted in Fig.JHJ along with the dependence of the gap on the chemical 
potential for a few values of g. This dependence makes the critical temperature sensitive to the presence of various 
impurities, complying with the reports of high-temperature superconductivity in doped graphite irP. 

These approximations are compared to the exact numerical solution of the self-consistency equation ( 37 1 , using 
the full Hamiltonian, Eq. ( 22 ) in Fig. [9] Moreover, Fig. 10 shows the position dependence of the gap away from the 
surface state for two different coupling strengths. For low values of the coupling it extends into the bulk only over 
a few interlayer distances due to a decay of the wave functions. Taking this into account we have chosen the model 
below Eq. (25), in which the order parameter is nonzero only on the outermost layers. 



D. Supercurrent 



Absence of dispersion in a flat band raises the questions of superconducting velocity and of the supercurrent: Can 
they be nonzero and, if they can, what is then the magnitude of the critical current? In this section we address the 
problem of supercurrent associated with the surface superconductivity in the flat-band multilayered rhombohedral 
graphene. Based on the model employed irP^ for description of the surface superconductivity the supercurrent was 
calculated irP^ as a response to a small gradient of the order parameter phase A = | A|e lkr using an approach similar 
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FIG. 9: Self-consistent surface gap function vs. coupling g. The blue line shows the results from the exact numerics in the case 
of TV = 20 layers and using /i — 0, the red line the numerical solution from Eq. ( |45[ ) at /i — and the magenta line at /J, = fj, opt 
plotted in Fig. [8] For large coupling, the surface gap tends towards the flat-band limit A oc g. For coupling g < Aira, the gap 
becomes exponentially suppressed, A oc exp(— AmaJg). The lower inset shows the ratio of A and the critical temperature T c as 
a function of the ratio g/a computed from Eq. ( |39| |. For low coupling it equals to the regular BCS value 1.764, whereas in the 
flat band regime g 3> 47ra it tends towards 3. 




Layer index 



FIG. 10: Profile of the gap function A for two different coupling strengths for TV = 20 layers. The lines are guides to the 
eye, as A is only defined in the discrete points marked with circles (g — 71) and squares (g = O.OI71). In the case g = 71 
the induced gap is already so large that high-momentum states with vfP ~ 71 contribute to superconductivity, and therefore 
superconductivity extends into the bulk of the material. 



to that o:P^ for supercurrent in a single layer of graphene. The supercurrent appears to be finite; the critical current 
is proportional to the superconducting zero-temperature gap, i.e., to the critical temperature, and to the radius of the 
flat band in the momentum space. Being produced by the surface superconductivity, the total current through the 
sample is independent of the sample thickness. Here we summarize the results obtained irP^Jfor the flat band regime. 
The current density along layer n is 

j„ = -ev F [ul{p)d-u n {p) + uj, (p)<ru n (p)] (1 - 2/ p ) . (46) 

P,9 



Here q labels different states for given p, while / p is the distribution function. The current operator in Eq. (46) 



couples the states at different sublattices and thus contains the overlaps of the wave functions localized at different 



surfaces. Let us consider a large but finite number of layers. Equation (17) tells us that the wave function for each 



sublattice decays away from the corresponding outer surface of the sample. As a result, the product of two wave 
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functions at different sublattices in Eq. (46) for the current, in addition to a decaying part, contains a contribution 
proportional to (pvf/ji) n = £ P /7i that does not depend on the distance from the surface. As is seen from Eq. (43) 
the energy £ p is of the order of the energy gap A at the surface. It is thus a characteristic energy associated with 
the superconducting coherence between the two surfaces. This coherence persists in the bulk even though the wave 
functions and the order parameter itself decay. This coherence produces a supercurrent that flows uniformly through 
the sample but with the density that is inversely proportional to the number of layers, see^. The full current through 
the sample is 

j _ 2ewAln(7i/A)k 
Trh 

where w is the width of the sample. The full current does not depend on the sample thickness Nd as expected 
for surface superconductivity. The critical current is determined by max(fc) ~ ^q -1 where the coherence length is 

ewAln(7i/A)p FB 



For nonzero fi we find in the same way as irP^ 

ewln(7i/A)k 



irti 



lAP 



(47) 



This equation holds for T <C |A|. 

The supercurrent is thus finite despite the absence of dispersion of the excitation spectrum. The critical current is 
proportional to the zero-temperature gap, i.e., to the superconducting critical temperature and to the size of the flat 
band in the momentum space. 



E. Effect of fluctuations 



In our analysis above we employ the mean-field approximation. The quality of this approximation is determined 
by the Ginzburg number which is a measure of the relative magnitude of order-parameter fluctuations. For usual 
3D superconductors the Ginzburg number is Gi ~ (T c /Ep) , it is very small due to a small ratio of the critical 
temperature to characteristic energy of electrons which is the Fermi energy. 

In the case of flat band surface superconductivity, both the critical temperature and the characteristic energy of 



electrons are of the same order. Indeed, according to results of Sec. IV B the characteristic energy £ p is of the order of 
A ~ T c . As a result, the Ginzburg number is of the order of unity. The magnitude of the order-parameter fluctuation 
Ai can be estimated as follows. The free energy density of fluctuations in the flat-band regime is 

F\ ~ A lPFB /h 2 . 

Since £ ~ fr/pFB one has the total free energy of fluctuations T\ ~ (?F\ ~ Ai. Comparing this energy with the 
thermal energy T we find Ai ~ T. The thermal fluctuations are thus of the order of the mean-field gap Ai ~ Ao for 
T ~ T c ; they freeze out only for T -C T c . Therefore, the mean-filed approach is not exact for the flat-band regime. It 
works well, however, for the regime when the quadratic spectrum dominates. 

If the superconductivity is dominated by the normal spectrum with a quadratic dispersion, the fluctuation free 
energy density for T not too close to T c is 

Fl ^^! =A 2 7071 
1 i i 



167r/i 2 Wp74 



where V2 — m* /2irh 2 is the 2D DOS and the effective mass m* is determined by Eq. (21). If the coherence length is 
£o = ^gAg" 1 where Ao is the mean-field gap, the full energy in an area 7r£g is 

2 Aj ^7o7i A2 747l 

Ag 16U|,74 Ag 7 
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Since Ji ~ T we find 

A i = Gi r , T ^o r , A o7o 
7471 7471 

If the quadratic dispersion dominates, the order parameter is A -C 7174/701 a s follows from Eq. ([2|. Therefore, 
the Ginzburg number Gi = e -1 ^ 2 <C 1, and the average fluctuation of the order parameter is small compared to its 
mean-field value. However, at the crossover point to the fiat band regime, Ao ~ 7174/70, and the fluctuation becomes 
of the same order as the mean-field value. 



F. Twinning boundary superconductivity 



The flat band states do not have to appear only at the surfaces of rhombohedral graphite. Similar flat bands 
can appear at twinning boundaries between different rhombohedrally stacked regions, or between a rhombohedral 
and Bernally stacked regiorP. The presence of such flat bands can be justified from the fact that at such twinning 
boundaries, the bulk topological charge Ni(p) discussed in^l changes sign or turns from a finite value to zero. Hence 
an asymptotically zero-energy state has to form at such boundaries for those momenta p in the p x ,p y plane for which 
Nx changes. 

We do not discuss the microscopic features of such twinning boundary flat bands here, but present only a numerically 
calculated gap function for such a system in a few example cases: (a) a stacking fault in rhombohedral stacking, where 
at one point in the sample the new layer stacked on top of the rhombohedrally stacked layers has its B-atom on top 
of the A atom of the previous layer, which on the other hand was on top of the B-atom of the previous layer (onset of 
Bernal stacking, but continued by rhombohedral stacking), (b) a few such Bernally stacked layers on top of each other, 
and (c) rhombohedral stacking fault in an otherwise Bernally stacked graphite. The corresponding nearest-neighbour 
Hamiltonians around the stacking fault region are of the form 

/H rhs (N 1 ,a+) \ 

71 (7- vpo ■ p 71 o + 

H a = 7i < T - v F cr ■ p 71 a- 

71°+ vpo ■ p 7i<7_ 

V H lhg (N 2 ,0 : / 

/H As (N 1 ,a+) \ 

7l<7_ VpO ■ p 7i<7 + 

7l<7_ Ufd • p 71 <7_ 

7i°+ v F a ■ p 71 er+ 

7l<7_ vpo ■ p 71 er_ 

V H lhg (N 2 ,o+)J 

(vpO-p 71 (7+ ^ 

7icr_ Vpo ■ p 71 er_ 

71 0+ vpo ■ p 71 o + 

//• H lhg (N,o + ) 

7lC+ VpO ■ p 71CT+ 

7ic_ vpo ■ p 7icr 

V 7i c+ v F<y ■ pj 

Here the dots denote zeros in the matrix and the matrices i? r h g (N, o) are the Hamiltonians for N layers of rhombo- 
hedrally stacked graphite with the coupling matrix o on the upper diagonal. 

The corresponding profiles of the superconducting order parameter are plotted in Fig. |11| This picture clearly 
shows how interface superconductivity shows up at the twinning boundaries. 



H h = 



V. SUMMARY 



The flat band with infinite DOS emerges in semi-metals with topologically protected nodal lines. This flat band 
promotes surface superconductivity with T c proportional to the pairing interaction strength and to the area of the flat 
band in the momentum space which is determined by the projection of the nodal line onto the surface. Topologically 
protected flat bands may also appear on interfaces, twin boundaries and grain boundaries in bulk 3D topological 
materials leading to an enhanced bulk T c . 
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FIG. 11: Profile of the gap function A for example systems showing twinning boundaries: (a) Two rhombohedrally stacked 
graphite (RHG) multilayers coupled by a stacking fault, (b) Two rhombohedrally stacked graphite multilayers coupled by a 
Bernally stacked multilayer consisting of 6 layers. In both (a) and (b) the upper rhombohedrally stacked part contains 20 
layers, and the lower part 10 layers, (c) Rhombohedrally stacked graphite multilayer sandwiched between two Bernally stacked 
multilayers. In all pictures, we have chosen g = O.I71 and disregarded the higher-order couplings 73 and 74. 



Rhombohedral graphite is a promising candidate for a system with topologically protected surface flat band at 
the Fermi energy. Earliej^ 10 ! 15 ! we have shown that, within the nearest-neighbour approximation, the rhombohedral 
graphite has a flat band for surface states, and these surface states support high-temperature superconductivity with 
the superconducting order parameter concentrated around the surfaces. The corresponding critical temperature is 
proportional to the pairing interaction strength and can be thus considerably higher than the usual exponentially 
small critical temperature in the bulk. This is in strong contrast to single-layer graphene, where the density of states 
vanishes at the Dirac point and therefore superconductivity quite generally requires strong dopinj22H241 j ( see a i so 
review^ and references therein) . 

However, next-nearest neighbour hoppings which are present in real rhombohedral graphite can break the exact 
topological protection and, therefore, the flat-band mechanism of superconductivity at sufficiently low values of the 
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coupling constant can be destroyed. Here we have studied the detailed effect of these higher-order interactions and 
demonstrated that instead of the flat-band scenario they can provide another mechanism of surface superconductivity 
which is of the BCS type but with a much larger coupling constant than the usual superconductivity in bulk. This 
large coupling constant comes from a large DOS associated with a heavy effective mass of surface quasiparticles that 
is clearly distinguishable on the background of the flat band which would exist without the higher-order interactions. 
For strong coupling, however, even this system crosses to the regime of flat band superconductivity. 

Indications towards surface superconductivit y hav e been seen in experiments on graphite in the form of a small 
Meissner effect and of a sharp drop in resistanc e! 16 * 17 !. The enhanced superconducting density has been also reported 
on twin boundaries in Ba^ei-zCo^AssP^. This year, these findings have been ratified and experimentalists have 
seen zero resistance in graphitic samples up to temperatures of 175 Kp^ and furthermore indications of even room- 
temperature superconductivity in specially prepared graphite samples^. These observations are compatible with 
surface or interface superconductivity described by our theory. However, they would require at least the presence of 
rhombohedrally stacked graphite regions embedded inside otherwise Bernally stacked regions of graphite. 

Besides the top and bottom surfaces of rhombohedral graphite, Bernally stacked graphite should have similar types 
of flat bands emerging on their lateral surfaces^^ZI. It is possible that these states also support high-temperature 
surface superconductivity. 

Our predictions provide a criterion for the parameters needed to obtain the highest critical temperature. They 
can be used for the search or for an artificial fabrication of layered and/or twinned systems with high- and even 
room-temperature superconductivity. 

We thank G. Volovik for helpful comments and the collaboration that initiated this project. We also acknowledge 
fruitful discussions with F. Mauri, A. Harju and M. Ijas. This work is supported in part by the Academy of Finland 
and its COE program 2012-2016, by the European Research Council (Grant No. 240362- Heattronics), and by the 
Program "Quantum Physics of Condensed Matter" of the Russian Academy of Sciences. 
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Note the difference between BdG holes that are time-reversed electron states and the valence band excitations (absence 
of electrons below the Fermi level), which are sometimes also referred as holes. By construction, the BdG equations are 
electron-hole symmetric. 



